For representation purposes, we limit ourselves to a rather small lattice size in the 3-dimensional case.
To extract the array of the matrix elements between eigenstates for the $\hat{n}_i = \hat{c}_i^\dagger\hat{c}_i$ operator at $i$-th site, we first need to extract the $i$-th row of the eigenvector arrays, since the eigenvectors are given as the columns of the said array as per scipy documentation. To obtain the array, we perform a tensor product operation. In our case, we select our site to be located at the center of the array.
$\langle \alpha|\hat{n}_i|\beta\rangle =\langle 0 |\sum_{j,k} \tilde{\alpha}^*_j\tilde{\beta}_k \hat{c}_j\hat{n}_i\hat{c}_k^\dagger|0\rangle = \tilde{\alpha}_i^*\tilde{\beta}_i$
To obtain the hopping matrix elements, we also extract the $i+m$ row of the eigenvector array. We then take the tensor product of the $i$-th and $(i+m)$-th row and add its conjugate (we also need to ensure all the conjugations in the tensor product are accounted for).
$\langle \alpha|\hat{h}_i^{(m)}|\beta\rangle = \langle \alpha|\hat{c}_i^\dagger\hat{c}_{i+m} + \hat{c}_{i+m}^\dagger\hat{c}_i|\beta\rangle = \tilde{\alpha}^*_{i}\tilde{\beta}_{i+m} + \tilde{\alpha}_{i+m}^*\tilde{\beta}_i$
Here, we provide two animations showing the eigenfunctions for the 2D model in the low and high disorder cases. Each frame in the corresponding animations displays one eigenfunction.